Massive Black Holes in Star Clusters. I. Equal-mass clusters 
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ABSTRACT 

In this paper we report results of coUisional A^-body simulations of the dynamical evolution 
of equal-mass star clusters containing a massive central black hole. Each cluster is composed of 
between 5,000 to 180,000 stars together with a central black hole which contains between 0.2% 
to 10% of the total cluster mass. 

We find that for large enough black hole masses, the central density follows a power-law 
distribution with slope p ^ t-^I-^s i^gide the radius of influence of the black hole, in agreement 
with predictions from earlier Fokkcr Planck and Monte Carlo models. The tidal disruption rate 
of stars is within a factor of two of that derived in previous studies. It seems impossible to grow 
an intermediate-mass black hole (IMBH) from a M < lOOM© progenitor in a globular cluster 
by the tidal disruption of stars, although M = 10^ Af© IMBHs can double their mass within a 
Hubble time in dense globular clusters. The same is true for the supermassive black hole at the 
centre of the Milky Way. 

Black holes in star clusters will feed mainly on stars tightly bound to them and the re- 
population of these stars causes the clusters to expand, reversing core-collapse without the need 
for dynamically active binaries. Close encounters of stars in the central cusp also lead to an 
increased mass loss rate in the form of high- velocity stars escaping from the cluster. A companion 
paper will extend these results to the multi-mass case. 

Subject headings: black hole physics — globular clusters — methods: N-body simulations — stellar dynamics 



1. Introduction 

Theoretical studies of the dynamics of massive 
black holes in dense stellar systems started in the 
1960s to explain the central activity and luminosi- 
ties of quasars. Since then, the dynamics of a mas- 
sive body in the centre of a stellar system has been 
the focus of a large number of theoretical stud- 
ies, starting with classic papers by Peebles (1972), 
Bahcall & Wolf (1976, 1977), Frank & Rees (1976) 
and Cohn & Kulsrud (1978). 

The problem is of great importance to astro- 
physics since the centres of the Milky Way and 
other nearby galaxies contain black holes of 10^ 
to IO'^Mq (Kormendy & Gebhardt 2001). In ad- 
dition, smaller sized black holes of a few thou- 
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sand solar masses might exist in globular clusters 
(Gerssen et al. 2002, 2003; Gebhardt et al. 2002; 
Portegies Zwart et al. 2004), although the evi- 
dence for them is still controversial (Baumgardt 
et al. 2003ab). A massive black hole in a galactic 
nucleus or a star cluster is a potential source of 
gravitational radiation due to its high mass and 
the fact that it will frequently undergo close en- 
counters with other stars and black holes if the 
density of the surrounding system of stars is large 
enough. It would therefore be a prime target for 
the forthcoming generation of ground and space- 
based gravitational wave detectors. 

Intermediate-mass black holes of several 100 to 
several 1000 M© could also be the explanation for 
the ultraluminous X-ray sources observed in ex- 
ternal galaxies (Portegies Zwart et al. 2004; DiS- 
tefano et al. 2004) and could provide the missing 
link between the stellar mass black holes formed 
as the end products of the stellar evolution of mas- 
sive stars and the 10^ — lO^Mp, sized black holes 



found in galactic centres (Ebisuzaki et al. 2001). 

Bahcall & Wolf (1976) showed that an equi- 
librium-flow solution for stars into the gravita- 
tional well around a black hole exists and pre- 
dicted that the stellar density will follow a power- 
law distribution p = r~" with exponent a — 1.75. 
While Peebles (1972) found a steeper slope, Monte 
Carlo simulations by Cohn & Kulsrud (1978) and 
Marchant & Shapiro (1980) confirmed the results 
of BahcaU & Wolf (1976). Due to the high stel- 
lar densities around the black hole, tidal disrup- 
tion of stars is important for the evolution of the 
system. Frank & Rees (1976) and Lightman & 
Shapiro (1977) found that stars on highly eccentric 
orbits dominate the consumption rate since stars 
drift faster in angular momentum space than in en- 
ergy space. Frank & Rees (1976) derived analytic 
formulae for the disruption rates which were later 
confirmed by Monte Carlo simulations (Marchant 
& Shapiro 1980; Duncan & Shapiro 1983). In the 
latter paper the authors also showed that black 
holes in globular clusters preferentially disrupt 
stars most tightly bound to the black hole while 
in galactic nuclei the disruption process should be 
dominated by stars not bound to the black hole. 
The gas lost from disrupted stars is either accreted 
onto the central black holes or lost in a stellar wind 
due to radiation drag, with the first process dom- 
inating for high enough black hole masses (David 
et al. 1987). 

Stellar collisions could also be an important 
process, although for the small mass black holes 
expected in globular clusters, Cohn & Kulsrud 
(1978) found that the collision rate of stars is 
about 30 times smaller than their consumption 
rate by the black hole. Similarly, Murphy et al. 
(1991) performed multi-mass Fokker-Planck cal- 
culations and found that in low-density galactic 
nuclei stellar disruptions happen more often than 
stellar collisions. 

Recently, Amaro-Seoane, Freitag & Spurzem 
(2004) followed the evolution of a system of equal- 
mass stars with a central black hole by means of 
an anisotropic gaseous model and found a strong 
black hole growth accompanied by the expansion 
of the cluster. An overall cluster expansion due to 
the inward drift and tidal disruption of stars was 
also predicted by Shapiro (1977). 

Although various aspects of the dynamical evo- 
lution of black holes in dense stellar systems have 



been studied in the literature, nobody has tried 
a self-consistent direct TV-body simulation of the 
growth of the central BH. Moreover, simulations 
with approximate methods such as Monte-Carlo 
or Fokker-Planck methods have been applied only 
to idealised systems. For example, with the ex- 
ception of Bahcall & Wolf (1977) and Freitag & 
Benz (2002), most simulations so far considered 
only single-mass systems, and ignored stellar evo- 
lution. Thus, though such studies are useful to 
gain physical understanding of the problem, they 
do not tell us much about the actual behaviour 
of star clusters with central black holes. Can an 
IMBH grow from a smaller mass seed black hole 
by accreting nearby stars? What will star clusters 
with an IMBH look like? Do they have cusps in 
surface luminosity? In the present and companion 
papers, we address these issues. Our focus wfll be 
mainly on the dynamics of star clusters contain- 
ing black holes. This is because direct TV-body 
simulations cannot be done for systems contain- 
ing more than a few times 10^ stars. In addition, 
simplifying assumptions like a fixed black hole at 
the cluster centre are most likely violated for stel- 
lar systems containing black holes of only a few 
hundred to a few thousand times the mass of a 
star. 

2. Description of the runs 

We simulated the evolution of star clusters con- 
taining between TV = 5,000 and N = 178,800 
stars using the coUisional Aarseth TV-body code 
NB0DY4 (Aarseth 1999) on the GRAPE6 boards 
of Tokyo University (Makino et al. 2003). AU clus- 
ters were treated as isolated and followed King 
profiles initially. At the start of the calculation, 
the massive black holes were at rest at the clus- 
ter centres. We set up the initial models so that 
the systems were in dynamical equilibrium and the 
density profiles were the same as the correspond- 
ing King models without the central black holes. 
To achieve this, we calculated the new distribution 
function from the King model density profile and 
the potential obtained by the original King poten- 
tial plus the BH potential, and generated positions 
and velocities of the stars using this new distribu- 
tion function. Note that it is impossible to setup 
an equilibrium model with isotropic velocity dis- 
persion, which has a flat core with finite central 
density around a black hole (Tremaine et al. 1994; 



Nakano & Makino 1999). Thus, our initial model 
is not in exact dynamical equilibrium, and the 
0.5% mass shell shows contraction of about 10% in 
case of a 5% BH mass. The effect is however much 
smaller than the initial contraction of the cluster 
with a central black hole in Fig. 5, which is caused 
by the development of an a = —1.75 cusp due to 
thermal evolution and can therefore be neglected. 

Two series of simulations were produced. In 
the first series we followed the evolution of equal- 
mass star clusters with central black holes. These 
simulations were designed to identify the relevant 
physical mechanisms and compare our results with 
theoretical estimates and results reported in the 
literature. In the second series of simulations we 
studied the dynamics of black holes in galactic 
globular clusters. We simulated multi-mass clus- 
ters containing between 16,384 < N < 131,072 
stars and a central black hole of Mbh = IOOOMq. 
The results of these runs will be reported in a com- 
panion paper (Baumgardt et al. 2004). 

In the present paper, stars were assumed tidally 
disrupted if their distance to the central black hole 
was smaller than a fixed disruption distance which 
was either rt = 10^^, 10^® or 10^^ in A^-body 
units. For solar type stars, the radius against tidal 
disruption by a 1000 M© IMBH is 2.3 • lO"'^ pc 
(Kochanek 1992, eq. 3.2). Since the half-mass 
radius of a star cluster is of order 1 in iV-body 
units and globular clusters have half-mass radii of 
several pc, the adopted tidal radii correspond ap- 
proximately to the tidal radius of an to = lAf© 
main-sequence star in a globular cluster. We as- 
sumed that a star was immediately disrupted if 
it entered the region with r < rt and was unaf- 
fected outside this area. The mass of a tidally dis- 
rupted star was added to the mass of the central 
black hole. To extrapolate the simulation results 
to larger N and compare with other results, we al- 
ways use fitting formulae which explicitly contain 
rt and scale them appropriately. 

So far our runs do not incorporate the effect of 
gravitational radiation, which should be unimpor- 
tant for the dynamical evolution of a star cluster as 
long as the stellar density around the black hole 
and its growth are dominated by main-sequence 
stars. 

In this paper, all clusters start from an ini- 
tial density profile given by a King model with 
Wq — 10.0 and are composed out of equal mass 



stars. Table 1 summarises other parameters. It 
first shows a number identifying the run and the 
number of cluster stars N. Shown next are the 
initial and final mass of the black hole divided 
by the mass of a single star, the tidal radius of 
the black hole, and the duration of the simula- 
tion. The last two quantities are given in A^- 
body units in which the constant of gravitation, 
total cluster mass and total energy are given by 
G ^ l,Mc ^ l,E ^ -0.25 respectively (Heg- 
gie & Matthieu 1986). The final columns contain 
a dimensionless constant describing the tidal dis- 
ruption rate of stars and the total number of tidal 
disruptions. We organised the runs into 4 groups. 
In the first set we varied the black hole mass and 
kept all other parameters constant. The next two 
groups contain runs where the number of cluster 
stars and the tidal radius was varied. The final 
group contains a few additional runs used in the 
paper. 

3. Theory 

Stars in the innermost cusp around the black 
hole, where the gravitational influence of the black 
hole is dominating, move on essentially keplerian 
orbits with slight perturbations when they en- 
counter other cluster stars. Bahcall & Wolf (1976) 
assumed that stars in the sphere of influence of 
the black hole follow an isotropic distribution in 
velocity space and are absorbed if their energy be- 
comes equal to the potential energy at the tidal 
radius. They then showed by Fokker-Planck calcu- 
lations that the stellar density distribution follows 
a power-law p{r) ~ r~" with a = 7/4 inside the 
sphere of influence of the black hole down to the 
radius where the tidal disruption of stars becomes 
important. The cusp profile will extend out to a 
radius where the self-gravity of the stellar system 
cannot be neglected any more. If the black hole 
mass is much smaller than the mass of the stars 
in the core, this happens when the velocity disper- 
sion in the cluster core becomes comparable to the 
circular velocity of stars in the field of the black 

hole: 

GMbh ,,. 

We found that 7 = 2 gives a good fit to the results 
of our simulations. The velocity dispersion of stars 
in the cluster core < v"^ > can be expressed in 



Table 1 


: Details of the 


performed A^ 


-body runs. 








No. 


N 


rl 


MBHi/m 


Mbh f/m 


^END 


ko 


ND^s 


1) 


80000 


10-7 


266.0 


827.0 


3000 


72.0 


561 


2) 


80000 


10-7 


800.0 


1388.0 


2000 


63.0 


588 


3) 


80000 


10-7 


2660.0 


3285.0 


2000 


63.1 


625 


4) 


80000 


10-7 


8000.0 


8749.0 


2000 


57.0 


749 


5) 


5000 


10-8 


800.0 


815.0 


2000 


44.1 


15 


6) 


10000 


10-8 


800.0 


823.0 


2000 


54.2 


23 


7) 


20000 


10-8 


800.0 


869.0 


2000 


65.0 


69 


8) 


80000 


10-8 


800.0 


997.0 


2000 


55.3 


197 


9) 


80000 


10-7 


800.0 


1388.0 


2000 


63.0 


588 


10) 


80000 


10-8 


800.0 


997.0 


2000 


55.3 


197 


11) 


80000 


10-9 


800.0 


851.0 


2000 


73.2 


51 


12) 


16000 


10-7 


2589.0 


2735.0 


2000 


75.6 


146 


13) 


20000 


10-7 


200.0 


338.0 


3000 


53.9 


138 


14) 


35700 


10-7 


1438.0 


1705.0 


2000 


63.0 


267 


15) 


65536 


10-8 


3276.0 


3513.0 


3000 


76.9 


237 


16) 


178800 


10-7 


461.0 


1368.0 


2000 


58.2 


907 



Notes: 

1. rt and Tend are given in A^-body units. 



terms of the core density and radius as: 



< C >= 



All GraUc r^ 



iVc 



= —Gmricr^ 



(2) 



Hence the influence radius is given by: 



3M 



BH 



STrmrirr^ 



15 r, 



Mbh 
m Nc 



(3) 



where Nc is the number of cluster stars and we 
assume that the core contains roughly 3% of all 
stars in the cluster. Assuming that the cusp profile 
goes over into a constant density core with density 
Uc at ri, the number of stars in the cusp can be 
estimated to be 



No 



usp 



An 



-n 

rin I — 



1.75 



r^dr 



250 ^»c 



(4) 



which serves to give an order of magnitude esti- 
mate. For central black holes containing less than 
a percent of the total cluster mass, the central 



cusp itself contains only of order 10% of the black 
hole mass in stars. Typical globular clusters would 
have black holes with masses of order 1000 Mq if 
they follow the relation found by Ferrarese & Mer- 
ritt (2000) and Gebhardt et al. (2000) for galaxies, 
making the detection of the central cusp in their 
density profile difficult even with HST (Drukier & 
Bailyn 2003). In clusters in which the black hole is 
more massive than the cluster core, an upper limit 
for Ti derives from the condition that the mass in 
stars inside r, should be smaller than the mass of 
the central black hole, i.e. M{< Vi) < Mbh- 

Frank & Rees (1976) found that rcrit, the radius 
at which tidal disruption of stars becomes impor- 
tant, is significantly larger than the tidal radius 
since relaxation lets stars drift faster in angular 
momentum space than in energy space. In order 
to account for this, they introduced the loss-cone, 
which is the area in angular momentum space con- 
taining all orbits with minimum distances smaller 
than the tidal radius rt of the black hole. At a 
given radius r, the opening angle of the loss-cone 
is given by 



.2 _ 2r, 



(5) 



for radii r < Vi. They then showed that the critical 
radius rcrit is approximately given by the radius 
where the time for stars to drift through the loss 
cone due to relaxation becomes longer than the 
crossing time. Inside Vcrit stars cannot drift in and 
out of the loss cone before falling into the black 
hole, so the loss-cone is empty. Cohn & Kulsrud 
(1978) and Marchant & Shapiro (1980) performed 
two-dimensional Monte Carlo calculations of the 
evolution of a star cluster with a central black hole 
which confirmed the loss-cone concept. They de- 
termined disruption rates of stars and showed that 
tidal disruption of stars flattens the density pro- 
file inside rcrit but has little influence outside this 
radius. 

Another process which can in principle be im- 
portant is the wandering of the black hole. A black 
hole in a stellar system is forced to move due to 
3 processes: Stars bound to the black hole force 
it to move around the common centre of gravity, 
stars escaping from the core cause a recoil motion 
of the black hole and the stellar core surrounding 
it, and passing unbound stars cause a brownian 
motion of the black hole in the centre of the clus- 
ter. Lin & Tremainc (1980) investigated the role of 
the different processes and concluded that passing 
unbound stars have the largest effect on the black 
hole. For a black hole in a constant density core, 
they estimated the wandering radius to be given 
by: 



is increased. 



Results 



^wand — U.y Tq 



Mbh 



(6) 



The wandering radius is difficult to determine in 
our simulations since it hardly exceeds the dis- 
tance of the innermost stars from the black hole 
and is therefore within the statistical uncertainty 
with which the position of the density centre can 
be determined. In addition, the innermost stars 
are so tightly bound to the black hole that they 
follow the motion of the black hole due to pass- 
ing unbound stars as long as the unbound stars 
pass at large enough distances. This further com- 
plicates the determination of the density centre. 
Chatterjee et al. (2002) gave a detailed discussion 
of black hole wandering and performed iV-body 
simulations which confirmed the validity of eq. 6. 
The wandering of the black hole will limit the for- 
mation of an a = 1.75 cusp to models with large 
enough black hole masses. It could also effect the 
capture rate of stars since the area of the loss cone 



4.1. Central density profile 

In order to determine the central density pro- 
files, we have overlaid 5 snapshots from the time 
the simulation was stopped for each cluster. All 
snapshots were centred on the black holes. Fig. 
1 depicts the final density profiles inside the half- 
mass radius for runs 1) to 4), which contain N ~ 
80, 000 stars and black holes with masses between 
300 < MBH/m < 8000. The influence radu de- 
fined by eq. 1 are marked by crosses. Also shown 
are the radii where the mass in stars becomes com- 
parable to the mass of the central black hole (solid 
circle) and the wandering radii of the black holes 
(open circle). For all cases studied, the critical 
radii are significantly smaller than the influence 
radii r^ and of the same order as the distance of 
the innermost stars from the black hole, so we can- 
not test the density profile inside Vcrit ■ The lack of 
stars inside Vcrit points to an efficient destruction 
at radii r < rcrit- 

For the black hole of a few hundred solar 
masses, an a = 1.75 cusp cannot be found with 
certainty since the number of stars inside r^ is too 
small. In addition, the wandering radius of the 
black hole is almost as large as r^, so the black 
hole wandering will also flatten the profile. Clus- 
ters with more massive black holes show a clear 
a = 1.75 cusp in good agreement with the predic- 
tion from Fokker-Planck and Monte Carlo models. 
The cusps extend up to radii r = r^ in all cases. 
As expected, r^ is significantly smaller than the 
radius where the mass in stars becomes compara- 
ble to the mass of the central black hole unless the 
black hole contains several percent of the cluster 
mass. If the black hole contains more than about 
5% of the total cluster mass, the a = 1.75 cusp 
goes all the way up to the radius where the mass 
in stars becomes equal to the mass of the central 
black hole. 

Fig. 2 depicts the central density profile for a 
number of small- A'^ runs. The overall behaviour is 
very similar to the high-A?^ runs, showing that our 
results do not depend on the particle number. Ex- 
tending our results to larger systems, we predict 
that in globular clusters only of order 50 to 100 
stars follow the central cusp profile if the mass of 




0-001 OOI 



Fig. 1. — Density profiles a,t T — 2000 A^-body units for clusters 1) to 4) from Table 1, which contain 
A^ = 80, 000 stars and black holes of varying masses. Solid circles mark the radii where the mass in stars 
becomes comparable to the mass of the central black hole, crosses mark r^ from eq. 1, open circles r^and 
and triangles show Vcrit- An a ~ 1.75 power-law cusp forms inside r^ in all models. For the model with the 
most massive black hole, the central cusp extends almost up to the radius where M(<r) ~ Mbh- 
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Fig. 2. — Same as Fig. 1 for clusters of lower mass. All clusters have a — 1.75 slopes inside the radii of 
influence of the black hole ri and flatter slopes outside this radius. 



the black hole is Mbh = IOOOMq. Since only a 
fraction of them will be bright enough to be easily 
detectable, it will be difhcult to find the central 
cusp in the luminosity profile of the cluster. In 
galactic nuclei with black hole masses in the range 
10*5 Mq < Mbh < lO^M© on the other hand, a 
considerable number of stars is in the a = 1.75 
cusp, so the detection of the black hole through 
observation of the central density or velocity pro- 
file should in principle be possible. 

4.2. Accretion rates 

We can estimate the rate at which stars are dis- 
rupted by the central black hole from the number 
of stars at radius Vcrit and the size of the loss cone, 
divided by the crossing time at Vcrit- 



.3/)2 



D 



n{r) 



TcR 



(7) 



Since in our simulations Vcrit is much smaller than 
ri, the density near the critical radius follows an 
a = 1.75 power-law distribution: 



n{r) = riQr 



(8) 



Following Frank & Rees (1976), the critical radius 
can be calculated to be: 



= 0.2 



rt^P 



BH 



rn^ no 



4/9 



(9) 



Hence we obtain for the disruption rate: 



D 



n m'J^ 



knVG 



p5/4 



r,^/%r/vw9 



M 



11/18 
BH 



(10) 



Fig. 3 shows the evolution of the disruption rate 
as a function of time for clusters with a range of 
black holes masses. In all runs the tidal disrup- 
tion rates decrease due to the cluster expansion 
which will be discussed in the next paragraph. The 
cluster expansion decreases the stellar density near 
rcrit and increases the crossing time, thereby de- 
creasing D. Solid lines show expected disruption 
rates according to eq. 10, calculated by determing 
no from the actual A'^-body data and the constants 
ko from a best-fit to the overall disruption rate. 



The time evolution of the disruption rates calcu- 
lated this way agrees very well with the one found 
in the iV-body simulations. In addition, the con- 
stants ko determined for the different simulations 
also agree reasonably well with each other and 
most are compatible with an average of kjj — 65 
(see Tabic 1). 

In order to calculate the collision rate in physi- 
cal units, we assume that the tidal radius of a star 
with radius i?, and mass m is given by: 



n 



1.3 «.(^) 



1/3 



(11) 



(Kochanek 1992, eq. 3.2). We then obtain for the 
tidal disruption rate: 



D 



0.00588 / R 



lOOMyrs V^ 



no 



pc- 



'V Mbh 



m 

WqJ V 1000^0 



(12) 



Although eq. 12 seems to indicate that D decreases 
with the black hole mass, this is not the case since 
the cusp density constant no also depends on the 
black hole mass and the central density Uc- It 
is however interesting to apply eq. 12 to cusps of 
observed systems. Genzel et al. (2003) for example 
derived a stellar density of p = 1.2 • IO^Mqpc"^ 
at a distance of r = 0.38 pc from the galactic 
centre. They found that the density inside this 
radius rises with a power-law with power a = 1.4, 
while it falls off with a ~ 2.0 outside this radius. 
Both values are not too far from the a = 1.75 
slope in our runs. Using their density, we find 
no = 2.21 • lO^pc^^-^^ if we assume an average 
stellar mass of <m>= IM©. From eqs. 9 and 12 
we then obtain a critical radius of rcrit — 1.2 pc 
and a total number of disruptions oi D ~ 30, 000 
per T — 100 Myrs for a central black hole mass of 
Mbh = 3 • 10^ Mq. Tidal disruption of stars could 
therefore play an important role for the growth of 
the galactic centre black hole. 

In order to apply eq. 12 to systems with a con- 
stant density core, we assume that the cusp den- 
sity goes over into a constant density core with 
density ric at r = 2ri. With the help of eq. 3 we 
then obtain 



D 



17504.9 
100 Myrs 



Rq 



4/9 



m 







-95/54 




T [NBODY] 



T [NBODY] 



Fig. 3. — Evolution of the tidal disruption rate of stars with time for 6 different clusters. Points with 
errorbars are results from A^-body simulations. Solid lines show a fit according to eq. 10 with the constants 
ko adjusted to match each case. The tidal disruption rate drops due to the expansion of the clusters. There 
is good agreement between the ko values for the different clusters. 
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lOOOMp 



61/27 

(13) 



The powers of the different factors are the same 
as the ones obtained by Frank & Rees for the case 
Tcrit < 'Ti (eq. 16a in their paper). Our disrup- 
tion rate is about a factor of 2 larger than theirs. 
From eq. 2, we can obtain an alternative expres- 
sion which uses the core velocity dispersion, which 
is easier to observe than the core radius: 
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This formula agrees with the one from Cohn & 
Kulsrud (1978) (eq. 66), who obtained nearly the 
same dependence of the disruption rate on the dif- 
ferent physical parameters. Our disruption rate is 
smaller by about 33%, which indicates a very good 
agreement given the errors involved by assuming 
that the central cusp goes over directly into a con- 
stant density core, which is a significant simplifi- 
cation of the real situation (Figs. 1 and 2). 

Typical parameters for the densest globular 
clusters are ric = 5 • lO^/pc"^ and Vc = 15 km/sec. 
While a Mbh = 1000 black hole should be able to 
double its mass within a few Gyrs from disrupted 
main-sequence stars in such a cluster, it does not 
seem possible to grow a 1000 Mq black hole from 
a 100 Mq progenitor within a Hubble time. A 
more detailed discussion must however also take 
the mass distribution of the stars into account, 
which will be the focus of a follow up paper. 

Fig. 4 shows the eccentricity distribution of 
stars which are tidally disrupted by the black hole 
and their semi major axis a in relation to the criti- 
cal radius. Stars that are disrupted by the central 
black hole move on very eccentric orbits on av- 
erage, with practically all of them having orbital 
eccentricities e > 0.999 on the final orbit prior to 
disruption, in good agreement with the idea that 
the drift in angular momentum space is more im- 
portant than the drift in energy space for the feed- 
ing of the black hole. The lower panel shows the 
semi major axis distribution. Most stars disrupted 
by the black hole have semi major axis a < Tcrit- 
The median a turns out to be about a = O.brcrit- 
Since the eccentricity is nearly unity, the apoboth- 
ron distance, i.e. the maximum black hole distance 
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Fig. 4. — Distribution of orbital eccentricities 
(above) and semi major axis (below) of stars being 
tidally disrupted by the central black hole for run 
2) which has N = 80,000 and Mbh = O-OlMci- 
All disrupted stars move on highly eccentric orbits, 
implying that drift in angular momentum space is 
the main process contributing to the tidal disrup- 
tion of stars. Stars being tidally disrupted by the 
black hole originate from approximately the criti- 
cal radius (bottom panel). 



on the last orbit before disruption is r^ = 2a, so 
Ta ~ Tcrit in good agreement with our numerical 
estimates for Tcrit- We do not find strong evi- 
dence for significant differences in the distribution 
of a/ Tcrit in the different runs, so we can average 
them. The dashed curve in the lower panel shows 
the distribution averaged over all runs, which is 
nearly identical to the solid curve for one particu- 
lar run. Since in our runs Tcrit « Ti all stars be- 
ing disrupted by the black hole are tightly bound 
to it. As these stars are constantly replenished by 
less bound stars from outside r^, energy conser- 
vation requires that the rest of the system gains 
energy and expands. This expansion will be the 
focus of the next section. 
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Fig. 5. — Evolution of lagrangian radii for two 
N = 65536 star clusters, one with a massive black 
hole and one without. The cluster without black 
hole goes into core-collapse at T = 1450 A^-body 
units and expands during the post-collapse phase 
due to energy generated by binaries formed in the 
core. The cluster with a black hole expands from 
the start since energy exchanges between stars in 
the cusp around the black hole provide the energy 
for the expansion. 
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Fig. 6. — Ratio of the lagrangian radii to the 
3% lagrangian radius for the run with a 1000 Mq 
black hole from Fig. 5 (solid lines). Dashed lines 
show average ratio near the end of the simula- 
tion. After an initial phase where the core is re- 
adjusting itself, the ratios of different radii become 
constant and the cluster expansion becomes self- 
similar. The black dots show the time when 5 ini- 
tial relaxation times have passed at the different 
radii. 



4.3. Cluster evolution 

The evolution of lagrangian radii of two clus- 
ters is depicted in Fig. 5. Shown is cluster 15), 
which contains N = 65536 stars and a black hole 
of Mbh = 0.05Mci and a second cluster with the 
same initial density profile and number of clus- 
ter stars but no massive black hole. The cluster 
without the black hole goes into core-collapse at 
T — 1460.0 iV-body units. Core-collapse is halted 
when binaries form in the core and interactions be- 
tween the binaries and passing cluster stars heat 
the cluster and lead to an overall expansion of the 
cluster. Since the cluster is isolated and we did not 
allow for stellar collisions, no characteristic length 
scale exists and the expansion is proportional to 



2/3 



(15) 



in the post-collapse phase (Giersz & Heggie 1994; 
Baumgardt et al. 2002). The cluster core under- 
goes gravothermal oscillations since the number 



of active binaries in the core is small: A look 
at the cluster data shows that most of the time 
there is only one binary in the core which pow- 
ers the expansion, in agreement with theoretical 
expectations (Goodman 1984). If this binary is 
expelled, the cluster centre re-collapses again and 
forms new binaries. In contrast, the cluster with 
a central black hole expands right from the start 
and without core-oscillations since the expansion 
is powered by energy exchanges of stars in the cen- 
tral cusp around the black hole, which remains 
in the cluster core due to its high mass. Only 
the innermost radii show a short collapse phase 
in the beginning, when the central cusp profile is 
created. Binaries cannot be responsible for the 
heating since their formation is suppressed by the 
high stellar velocities in the cusp and the strong 
gravitational field of the black hole which disrupts 
binaries. Consequently, no stable systems formed 
in our runs. Fig. 5 confirms that a central black 
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hole can act as a heat source similar to the effect 
of binaries in post-collapse clusters. 

Fig. 6 shows the evolution of lagrangian radii 
for cluster 15), which has A^ = 65536 stars and a 
black hole mass of O.ObMci- AH radii are divided 
by the 3% lagrangian radius. If we calculate the 
local relaxation time at different lagrangian radii 
from Spitzer (1987): 



tr = 0.065 



nm^lnA 



(16) 



with A = O.lliV (Giersz & Heggie 1994), it can 
be seen that the ratio of the different radii to the 
3% radius becomes constant after about 5 local 
relaxation times have passed. The expansion of 
the cluster becomes therefore self-similar beyond 
this time. Since the relaxation time increases with 
radius, the equilibrium profile is established first 
for radii nearest to the black hole and then forms 
at larger radii. We obtain similar results for the 
other clusters. Most globular clusters have cen- 
tral relaxation times much less than a Hubble time 
(Spitzerl987, Fig 1.3), so we expect that they have 
reached equilibrium profiles in their centres if they 
contain massive black holes. 

Closely related to the expansion of a cluster is 
the escape of stars. The energy distribution of 
stars escaping from the two clusters of Fig. 5 is 
depicted in Fig. 7. We measure the energies of 
escaping stars after they have left the cluster and 
divide the energies by the average kinetic energy 
of all stars still bound to the cluster at the time 
the escape event happens. The energy distribution 
of escapers in case of no black hole shows two dis- 
tinct maxima, corresponding to escapers created 
by a slow diffusion process in the outer parts of 
the cluster and escapers created by three-body en- 
counters in the cluster core. We can fit the whole 
escaper distribution by the sum of two gaussians 
with maxima at < log E/EKin >= 0.1 and 10.0, 
similar to what Baumgardt et al. (2002) found for 
isolated clusters with lower N. 

Although the structure of both clusters is not 
exactly the same and a direct comparison therefore 
difficult, it can be seen that the number of escapers 
is increased if a black hole is present (note differ- 
ent scale on y-axis in lower panel). At T=1450, 
at which time the cluster without IMBH reached 
core collapse, the cluster with a black hole has 
created about 6.5 times as many escapers as the 
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Fig. 7. — Distribution of escape energies for the 
clusters of Fig. 5. For a cluster without central 
black hole, two distinct peaks are visible, corre- 
sponding to stars escaping due to distant two- 
body encounters and close three-body encounters 
involving binaries. In the cluster with a black hole, 
most stars escape due to close encounters in the 
cusp around the black hole. Solid lines are gaus- 
sian fits to the energy distributions of different es- 
cape mechanisms (see text). Dashed lines show 
the sums, which fit the A''-body results very well. 

cluster without black hole. In the post-collapse 
phase, the number of escapers per time is approx- 
imately equal in the two runs. Since the process 
for generating escapers is different, the energy dis- 
tribution of escapers in the two runs is also differ- 
ent especially at the high-E end. Escapers created 
by distant encounters outside the core should still 
be present in the black hole case. If we subtract 
the distribution found for stars escaping due to 
distant encounters in case of no black hole from 
the black hole case, we can fit the remaining dis- 
tribution by a gaussian distribution with mean 
< logE/Exin >~ 1.0. These are escapers cre- 
ated in the high-density cusp around the black 
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hole. From Fig. 7 we can estimate that about 200 
stars escape by distant encounters and about 420 
by close encounters in the cusp. During the same 
time 211 stars are disrupted by the black hole, so 
stars in the cusp have a larger chance of escape 
than disruption. Nevertheless, the energy carried 
away by the escapers is small compared to the en- 
ergy created by the tidal disruption of stars, since 
we find that about 10 times as much energy is cre- 
ated than taken away by escapers. The energy 
created by the tidal disruption of stars therefore 
enters mostly the cluster expansion. We obtain 
this result in all simulated clusters. 

Mathematically, when a star is accreted to the 
central BH, the binding energy of the star cluster 
is reduced by the amount same as the binding en- 
ergy of the star and the IMBH. In other words, 
the cluster is heated up. Since almost all stars 
accreted to BH are strongly bound to IMBH, ac- 
cretion almost always resulted in heating. From 
a physical point of view, accreted stars increased 
their binding energy by giving their kinetic ener- 
gies to field stars, thus heating the cluster. How 
much heat the accreted star gave is simply deter- 
mined by its binding energy at the moment of the 
accretion. The effective energy generation by the 
stars disrupted by the black hole can be estimated 
from the disruption rate multiplied by the typi- 
cal energy of a star disrupted by the black hole. 
Since Vcrit » ft, stars disrupted by the black 
hole are on nearly radial orbits, so the energy of 
a star is equal to the potential energy at rcrit' 
E^ = G MBHiTi/rcrit- We therefore obtain 



E ^ 



D-E, 



M 



1/2 
BH 



(17) 



Interestingly, there is no dependence of the en- 
ergy generation rate on the tidal radius r^. For 
small mass black holes, most of the potential en- 
ergy comes from the gravitational interaction be- 
tween the stars themselves. The energy needed to 
expand the cluster is therefore given by 



E^k^^ 



Th 



(18) 



with r^ being the cluster's half-mass radius and k 
a constant of order unity. Combining both equa- 
tions and assuming that the cusp profile goes over 
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Fig. 8. — Anisotropy as a function of radius for 
two clusters with A^ = 80, 000 stars. The dotted 
vertical line marks the radius of influence of the 
black hole. The cluster with the larger black hole 
has a slight tangential anisotropy in the centre due 
to the preferential disruption of stars on radial or- 
bits. 

into a constant density core at r^, we obtain 





(19) 


With ko = 65, Mci = 1.0 and rh/rc w 


10 this 


gives 




rT-r%' ^6-10'' Ml^j^mt 


(20) 



The time dependence is the same as in the self- 
similar case, r ^ t^", since the energy generation 
rate is independent of the tidal radius, so no char- 
acteristic length scale exists. 

Fig. 8 depicts the velocity anisotropy profile for 
runs 2) and 4), which have N — 80,000 clusters 
and black holes containing 1% and 10% of the total 
cluster mass. The anisotropics were defined as 



/?=! 






(21) 
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where Vr and vt are the radial and tangential ve- 
locities of each star and the sum runs over all stars 
in a radial bin. No significant degree of anisotropy 
can be detected in the cluster with the small black 
hole inside the influence radius of the black hole. 
The cluster with the more massive black hole 
has a slightly tangentially anisotropic velocity dis- 
persion, similar to what Cohn & Kulsrud (1978) 
found in their Fokker-Planck calculations for stars 
inside the influence radius of the black hole. The 
reason is the preferential disruption of stars on ra- 
dial orbits and the fact that more stars are dis- 
rupted in a cluster with a more massive black hole. 
Amaro-Seoane, Freitag & Spurzem (2004) (their 
Fig. 8) obtained a radially anisotropic profile be- 
tween Tcrit < r < Ti, which could be due to their 
very large half-mass radius at the end of the run, 
which decreases the merging rate of stars. In any 
case, the small amount of anisotropy together with 
the small number of stars inside r = Vi makes an 
observation of this effect for globular clusters al- 
most impossible, although it might be detectable 
in galactic nuclei. For isolated clusters, cluster ha- 
los are build up from stars scattered out of the cen- 
tre which move on very radial orbits. Therefore, 
both clusters have radially anisotropic profiles out- 
side the half-mass radius. For any realistic cluster, 
the tidal field would cause an isotropisation of the 
stellar orbits. 

5. Conclusions 

We report the first results of self-consistent 
A'^-body simulations of star clusters composed of 
equal-mass stars and a central massive black hole. 
We find that in clusters with a massive central 
black hole a p oc r^^-"^^ power-law cusp forms 
inside the sphere of influence of the black hole, 
in good agreement with predictions from Fokker- 
Planck and Monte Carlo simulations. In star clus- 
ters where the black hole mass is less than a few 
percent of the total cluster mass, the cusp contains 
only a fraction of the black hole mass in stars. 
The minimum black hole mass needed to form an 
a = 1.75 cusp is several lOOM© for a typical glob- 
ular cluster. Otherwise the cusp contains too few 
stars to be significant. For a black hole mass less 
than a few percent of the total cluster mass, the 
density profile is shallower than a = 1.75 outside 
the radius of influence r^ of the black hole. For 
more massive black holes, the a = 1.75 cusp ex- 



tends all the way through the core. 

The cusp proflle forms from the inside out and 
it takes about 5 local relaxation times until it is 
established at a given radius. Since central relax- 
ation times of globular clusters are of order 10^ to 
10^ years, globular clusters should have a = 1.75 
power-law cusps in their centres and will evolve 
more or less along a sequence of equilibrium pro- 
files if they contain massive black holes. Inside 
the radius of influence of the black hole, the veloc- 
ity profile is slightly tangentially anisotropic due 
to the tidal disruption of stars on radial orbits. 
Since the magnitude of this effect is small, it is 
unlikely that it will be detectable except for star 
clusters with very massive black holes. 

Our simulations confirm the merging rates 
found in Fokker-Planck simulations and from an- 
alytic estimates. Tidal disruption of stars could 
play an important role for the current growth of 
the super-massive black hole at the galactic centre 
and the growth of intermediate-mass black holes 
in dense star clusters. However, the formation of 
an intermediate-mass black hole out of a stellar 
mass one in a globular cluster by tidal disruption 
of stars alone seems impossible. A massive black 
hole in a star cluster merges mainly with tightly 
bound stars from its direct vicinity. These stars 
are constantly replaced by less bound cluster stars 
which drift inward due to relaxation. This causes 
an overall expansion of the cluster. Black holes 
can therefore halt the core collapse of globular 
clusters, similar to the effect of binaries in post- 
collapse clusters. As in the case of an expansion 
driven by a central population of binaries, the 
expansion is self-similar, i.e. r ^ t^'^ . 
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